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Nomenclature 
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elementary  charge 

Us 

=  ion  acoustic  velocity 

k 
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Boltzmann  constant 

Sw 

=  secondary  electron 

ne 

= 

electron  number  density 

emission  yield 

nt 
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ion  number  density 

Ad 

=  Debye  length 

n 
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non-dimensional  ion  number  density 

(p 

=  potential 

me 
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electron  mass 

1 

=  non-dimensional  potential 

nij 

= 

ion  mass 

Te 

= 

electron  temperature 

= 

ion  velocity 

U 
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non-dimensional  ion  velocity 

u 

= 

axial  component  of  the  non-dimensional  ion  velocity 

V 

= 

radial  component  of  the  non-dimensional  ion  velocity 
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I.  Introduction 

Recent  experimental  studies  have  revealed  the  presence  of  high  velocity  ions  at  high  angles  off  the  centerline  in 
Hall  thruster  plumes  [  1  ] ,  [2] .  These  ions  have  kinetic  energies  that  are  higher  than  can  be  explained  by  collisional 
processes  such  as  charge  exchange  and  elastic  scattering,  and  are  thought  to  be  ions  that  have  been  accelerated 
through  nearly  the  full  acceleration  potential  of  the  thruster  [3].  From  a  thruster  integration  standpoint,  such  ions  are 
of  definite  concern  because  they  can  cause  damage  to  surfaces  with  which  they  collide.  Therefore,  it  is  important  to 
understand  the  underlying  mechanism  which  generates  these  ions. 

We  hypothesize  that  the  radial  electric  field  necessary  to  bend  ion  trajectories  towards  high  angles  could  be 
provided  by  the  sheath  at  the  corner  of  the  Hall  thruster  exit,  as  illustrated  in  Figure  1.  In  a  quasi-neutral  plasma,  a 
region  of  non-neutrality  will  develop  naturally  at  a  wall  due  to  the  imbalance  between  the  thermal  electron  and  ion 
fluxes  towards  that  surface.  This  non-neutral  region,  the  sheath,  has  its  own  internal  electric  field  that  is  directed 
towards  the  wall.  So  the  question  becomes:  is  the  electric  field  provided  by  the  sheath  sufficient  to  accelerate  ions 
with  a  certain  initial  axial  velocity  to  a  certain  angle?  To  answer  this  question,  the  thickness  and  potential  profile  of 
the  sheath  must  be  determined.  Additionally,  to 
apply  this  problem  to  a  Hall  thruster,  the  fact  that 
the  wall  consists  of  a  dielectric  material  must  also 
be  taken  into  account. 

Previous  research  has  addressed  the  problem  of 
sheath  development  in  a  flowing  plasma.  Hong  and 
Emmert,  for  example,  have  investigated  sheath 
behavior  in  the  wake  of  a  metal  target  [4].  While 
useful  from  a  problem  formulation  standpoint  (the 
same  equations  and  numerical  solution  methods  can 
be  applied  to  dielectric  corner  problem),  the 
geometries  examined  in  these  cases  were  not  similar 
to  that  of  current  problem,  nor  were  the  effects  of  a 
dielectric  material  considered. 

Studies  conducted  by  Ahedo  have  treated  the 
effects  of  plasma  flowing  past  annular  dielectric 
walls,  as  is  the  case  in  a  Hall  thruster.  In  Ahedo's  work,  a  model  was  first  developed  to  describe  the  pre -sheath  (the 
quasi-neutral  region  outside  the  sheath)  [5].  This  model  was  then  linked  to  equations  describing  the  non-neutral 
sheath,  so  that  ultimately  the  entire  domain,  from  one  wall  to  the  other,  was  simulated.  Secondary  electron  emission 
at  the  walls  was  also  taken  into  account  [6].  However,  this  approach  cannot  be  directly  applied  to  the  problem  of 
radial  acceleration  within  the  sheath,  because  it  assumed  a  "zero-Debye  length  limit,"  i.e.  the  sheath  was  assumed  to 
have  zero  thickness.  Also,  the  presence  of  the  corner  at  the  exit  adds  a  further  complication  to  the  direct  application 
of  this  method  to  the  problem  at  hand. 

To  determine  whether  the  sheath  at  a  dielectric  corner  has  a  significant  effect  on  ion  trajectories,  a  2D  framework 
has  been  developed  to  solve  the  standard  sheath  equations  while  including  the  effects  of  secondary  electron 
emission.  Ultimately,  the  goal  is  to  apply  this  2D  model  to  the  entire  region  around  a  dielectric  corner  in  a  flowing 
plasma.  This  paper  outlines  the  steps  carried  out  so  far  in  this  effort,  including  the  implementation  and  validation  of 
the  equation  solver,  as  well  as  the  application  of  the  model  in  a  simplified  "quasi"  2D  fashion  to  obtain  first 
estimates  of  sheath  thickness,  potential  profile,  and  radial  electric  field. 

II.  Model  Formulation 

To  model  the  electric  field  within  the  sheath,  a  2D  code  was  developed.  This  code  solves  continuity, 
momentum,  and  Poisson’s  equations  iteratively  to  determine  density,  velocity,  and  potential  profiles  within  the 
sheath.  The  direction  perpendicular  to  the  wall  normal  vector  will  be  called  the  "axial"  or  "z"  direction,  while  the 
direction  parallel  to  the  wall  normal  vector  will  be  called  the  "radial"  or  "r"  direction  (to  correspond  to  Hall  thruster 
convention). 

A.  Sheath  Equations 

To  model  the  sheath,  the  familiar  hydrodynamic  sheath  equations  were  applied  [4].  These  equations  consist  of 
ion  continuity  and  ion  momentum  conservation,  as  well  as  Poisson's  equation.  Electron  density  is  modeled  using  the 
Boltzmann  relation.  The  equations,  which  hold  for  a  collisionless  plasma,  are  as  follows: 
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Figure  1  -  Possible  bending  of  ion  trajectories  due  to  the 
sheath.  Ions  with  initially  axial  trajectories  may  be  deflected 
by  the  sheath's  internal  electric  field. 
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To  simplify  the  calculation,  one  can  combine  equations  (3)  and  (4)  and  introduce  the  following  non-dimensional 
parameters: 
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(the  Debye  length  and  the  ion  acoustic  speed,  respectively).  Substituting 
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B.  Numerical  Method 

The  problem  was  divided  into  two  separate  parts:  a  fluid  solver  and  a  potential  solver.  The  fluid  solver 
determines  the  values  of  ion  density  and  velocity  from  the  continuity  and  momentum  equations,  given  a  potential 
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field.  The  potential  solver  determines  the  potential  field  from  Poisson’s  equation  based  on  the  ion  density  and  the 
potential  from  the  previous  time  step. 

A  finite  volume  approach  was  used  to  calculate  the  spatial  derivatives  of  the  fluid  equations  on  a  rectangular  grid 
with  uniform  spacing  in  the  axial  and  radial  directions.  First  order  upwind  values  were  used  to  compute  velocities  at 
the  boundaries  of  each  cell,  while  density  was  calculated  at  the  cell  centers.  Time  derivatives  were  calculated  using 
the  Beam-Warming  algorithm  to  ensure  stability  [4], [7].  Finite  volumes  were  also  applied  to  the  spatial  derivatives 
of  Poisson’s  equation,  and  the  potential  field  was  solved  by  applying  an  iterative  Gauss-Seidel  scheme  [8]. 

C.  Model  Verification 

The  model  was  validated  by  comparing  its  result  to  those  obtained  using  the  standard  ID  sheath  equations  [9]. 
In  ID,  the  non-dimensional  sheath  equations  are: 


nU  =  1 

(12) 

U2  +  2%  =  1 

(13) 

dr  ^ 

—  =  -C-exp^_  (14) 

Op 


Where  p  —  r!/ iD  .  By  using  (12)  and  (13)  to  solve  for  the  ion  density,  one  can  derive  the  following: 


d2  r  ^ 

—  =  -<-2j-exp^2  (15> 

dp 


Equation  (15)  is  an  ODE  that  can  be  solved  using  standard  methods  given  conditions  for^  and  d%/dp  at  the 
outer  boundary,  and  either  the  sheath  thickness  (if  known)  or  the  value  of  %  at  the  wall  boundary.  Figure  3  and 

Figure  4  shows  a  comparison  of  the  ID  solution  and  the  results  of  the  2D  solver.  In  each  case,  % ,  d%j dp ,  and 

,  were  used  to  determine  the  potential  profile  within  the  sheath  and  the  sheath  thickness.  The  comparisons 

shown  in  Figure  3  and  Figure  4  suggest  that  the  2D  code  solves  the  sheath  equations  with  an  acceptable  degree  of 
accuracy,  and  thus  can  be  applied  to  the  problem  of  the  sheath  at  a  dielectric  corner. 


D.  “Quasi”-2D  Application 

The  intent  is  ultimately  to  apply  the  2D  solver  to 
a  domain  that  encompasses  the  entire  area 
surrounding  the  corner,  as  shown  in  Figure  2. 
However,  to  carry  out  the  calculation  in  this  region, 
a  major  obstacle  must  be  overcome.  It  is  important 
to  note  that  as  the  value  of  the  potential  on  the  outer 
boundary  changes,  so  does  the  thickness  of  the 
sheath.  Therefore,  because  the  2D  solver  is  only 
valid  in  the  region  inside  the  sheath,  in  order  to 
implement  the  solver  in  its  current  form,  the  width 
of  the  simulation  domain  must  be  allowed  to  vary. 
Alternatively,  rather  than  having  a  variable-width 
domain,  one  could  attempt  to  solve  the  sheath  as 
well  as  the  quasi-neutral  region  adjacent  to  it  (the 
"pre-sheath").  This  has  been  done,  for  instance,  by 
Ahedo  [5], [6].  Additionally,  techniques  for  joining 
the  two  regions  have  been  developed,  such  as 
asymptotic  matching  [10]  and  patching  [11]. 


solution  domain 


quasi-neutral  region  ^ 

- _ 

rs 

■ - ~~~shea 

th 

dielectric  corner 

Figure  2  -  Proposed  simulation  domain.  If  the  entire  pre¬ 
sheath/sheath  solution  is  obtained,  the  solution  domain  must 
be  split  up  into  two  regions.  If  just  the  sheath  solution  is 
obtained,  the  width  of  the  domain  must  vary. 
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Normalized  Potential  from  2D  solver,  £=  et()/kTe 


0  0.2  0.4  0.6  0.8  1 


Comparison  of  1 D  and  2D  solutions 


i_ i_ i_ i_ i_ i_ i 

0  0.5  1  1.5  2  2.5  3  3.5 

P  =  r^D 

Figure  3  -  ID  solution  vs.  2D  solver  results,  Example  1.  In  this  case  %  =  0,  and  d%  / dp  =  -0.1  at 
the  outer  boundary,  and  %w  =-0.4. 


Comparison  of  ID  and  2D  solutions 


Figure  4  -  ID  solution  vs.  2D  solver  results,  Example  2.  In  this  case  %  =  0,  and  d%/ dp  =  -0.1  at 

the  outer  boundary,  and  =-0.8  5 
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Rather  than  start  by  attempting  to  implement  a  variable-width  domain  or  solve  the  entire  pre-sheath/sheath 
problem,  an  intermediate  step  was  taken.  In  this  step,  the  2D  solver  was  applied  to  small  cells  along  the  corner.  The 
cells  were  sized  such  that  it  could  be  assumed  that  the  potential,  electric  field,  density,  and  velocity  were  constant 
along  the  outer  boundary  of  the  cell,  as  illustrated  in  Figure  5.  Then  for  each  cell,  the  2D  solver  was  used  to  find  the 
variation  of  the  potential  within  the  sheath,  as  well  as  the  sheath  thickness.  In  this  fashion,  the  variation  of  the 
sheath  thickness  and  the  electric  field  in  the  axial  direction  could  be  estimated. 


Figure  5  -  Applying  the  2D  solver  in  a  "quasi"  2D  fashion.  The  solver  was  applied  to  cells  over 
which  the  outer  boundary  values  could  be  assumed  constant. 


The  boundary  conditions  for  each  individual  cell  are  outlined  in  Figure  5.  At  the  outer  boundary,  the  value  of 
%  ,  d%  I  dp  ,  n  and  U  were  set  at  fixed  values.  Note  that  while  %  and  d%  I  dp  differ  depending  on  the  cell's 

axial  location,  n  —  1  and  U  —  1  at  the  edge  of  the  sheath  (since  nf  =  n0  and  w,  =  us  here).  The  sheath  thickness  was 
then  was  varied  until  the  value  of  %  at  the  wall  satisfied  the  condition: 


<-<V 


(16) 


Where  8W  =  8W  ^  ^  is  the  secondary  electron  emission  yield,  an  empirical  relationship  that  describes  the 

proportion  of  secondary  electrons  generated  when  a  primary  electron  collides  with  a  particular  material.  For  Hall 
thrusters,  this  material  is  typically  Boron  Nitride  or  Borosil.  However,  it  should  be  noted  that  for  high  values  of  Te , 
as  are  typically  found  near  the  exit  plane  of  a  Hall  thruster,  the  sheath  reaches  its  "charge  saturation  limit". 

Therefore,  for  the  purposes  of  the  "first  estimate"  detailed  in  this  paper  it  will  be  assumed  that  %w  ~  -1*018  (the 
charge  saturation  value). 

When  applicable,  non-dimensional  boundary  values  were  taken  from  an  HPHall  model  of  an  SPT-70  Hall 
thruster.  HPHall  is  a  2D  hybrid-PIC  model  in  which  the  electrons  are  treated  as  a  fluid  and  the  ions  are  treated  as 
particles  [12].  The  particular  version  of  the  code  that  was  used  was  one  that  has  been  under  development  at  JPL  (see 
[13]).  The  SPT-70,  a  Russian-designed  thruster,  was  selected  because  of  the  wide  number  of  experimental  studies 
that  have  been  conducted  on  it,  as  well  as  the  fact  that  it  has  been  used  to  vet  Hall  thruster  codes  in  the  past  [12]. 
Conditions  on  %  and  d%!  dp  for  the  SPT-70  at  the  outer  boundary  of  the  solution  domain  are  shown  in  Figure  6. 
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Non-dimensional  Potential  (%)  at  Outer  Boundary 


Non-dimensional  Radial  Electric  Field  (d%/dp)  at  Outer  Boundary 


C  =  z/XD 


Figure  6  -  Boundary  conditions  from  the  HPHall  SPT-70 
solution.  Note  that  £  =  0  corresponds  to  the  thruster  exit  plane. 

III.  Results  and  Discussion 

Applying  the  boundary  conditions  shown  in  Figure  6  according  to  the  method  outlined  in  the  preceding  section 
produce  the  results  shown  in  Figure  7.  These  results  clearly  illustrate  that  the  sheath's  thickness  increases  as  the 
distance  from  the  exit  plane  is  decreased.  Additionally,  the  potential  drop  across  the  sheath  decreases  as  one  moves 
closer  to  the  exit,  resulting  in  a  weaker  radial  electric  field  near  the  exit  plane.  Near  the  exit,  the  strength  of  the 
radial  electric  field  is  about  an  order  of  magnitude  smaller  than  it  is  at  the  upstream  boundary. 

These  results  suggest  that  an  axially-moving  ion  may  start  outside  the  sheath,  but  as  it  moves  towards  the  exit 
plane  it  can  enter  the  sheath  and  be  accelerated  by  the  radial  electric  field.  However,  these  results  also  indicate  that 
the  radial  electric  field  strength  decreases  significantly  towards  the  exit  plane.  This  means  there  is  a  trade-off 
between  where  the  ion  enters  the  sheath  and  the  amount  of  radial  acceleration  that  it  experiences.  The  angle  at 
which  the  ion  enters  the  sheath  also  influences  its  subsequent  trajectory.  Therefore,  to  understand  the  true  effect  of 
the  sheath  on  ion  trajectories,  several  "test"  particles  were  introduced  at  the  sheath  boundary  and  their  subsequent 
paths  tracked. 

Figure  8  to  Figure  10  shows  the  results  of  ion  tracking  within  the  sheath.  Each  plot  represents  a  different  initial 
axial  velocity.  The  sheath  edge  is  plotted  in  red,  while  the  particle  paths  originating  from  several  different  starting 
points  are  shown  in  blue.  Table  1  to  Table  3  catalog  starting  and  ending  values  for  the  trajectories  shown  in  these 
figures. 
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Figure  8  show  the  results  for  an  initial  axial  velocity  of  U0  =  5  (note  that  V0 ,  the  initial  radial  velocity,  is  set  at  -1 
to  correspond  to  the  Bohm  velocity  into  the  sheath).  As  can  be  seen  from  this  figure,  two  of  the  trajectories  intersect 
the  wall,  while  the  others  pass  the  exit  plane  without  colliding  with  the  corner.  Of  the  two  trajectories  that  do  not  hit 
the  wall,  trajectory  (3)  has  a  change  in  angle  of  about  1 1°,  while  trajectory  (4)  has  a  change  in  angle  of  about  2°. 

Looking  at  Figure  9,  one  can  see  that  for  U0  -  10,  only  trajectory  (1)  intersects  the  wall,  while  the  other 
trajectories  pass  through  the  exit  plane.  In  this  case,  as  shown  by  Table  2,  the  angle  change  for  the  unobstructed 
trajectories  ranges  from  roughly  0.5  (trajectory  (4))  to  5°  (trajectory  (2)).  In  the  case  of  U0  =20  (Figure  10  and  Table 
3,  all  trajectories  pass  through  the  exit  plane,  but  the  angle  change  ranges  from  about  0  to  2°. 

These  data  show  that  the  slower  the  initial  axial  speed,  the  larger  the  change  in  angle  that  can  be  produced. 
However,  if  the  axial  speed  is  too  low,  or  the  ion  enters  the  sheath  too  far  upstream  of  the  exit  plane,  then  the  ion 
will  collide  with  the  wall  rather  than  leave  the  thruster.  Note  that  the  for  an  SPT-70  thruster  operating  at  nominal 
conditions,  U0  is  approximately  equal  to  5.  Therefore,  one  would  expect  to  see  a  maximum  angle  change  on  the 
order  of  10°  due  to  the  sheath  in  an  SPT-70.  This  is  substantially  less  than  the  angle  change  required  to  produce  the 
high  velocity,  high  angle  population  seen  in  experiments.  To  reproduce  the  experimental  results,  one  would  have  to 
see  a  change  in  angle  of  60  to  80°  [1],[2].  This  suggests  that  another  mechanism  may  be  responsible  for  producing 
high  velocity,  high  angle  ions. 


Non-dimensional  Potential,  ecfATe 


Non-dimensional  Radial  Electric  Field,  3^/Sp 


Figure  7  -  Non-dimensional  potential  and  non-dimensional  electric  field,  plotted  as  a  function 
of  axial  and  radial  location.  The  white  area  represents  the  quasi-neutral  region  outside  the  sheath. 
Note  that  £  =  0  corresnonds  to  the  thruster  exit  nlane. 
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Figure  8  -  Ion  trajectories  originating  from  four  different  starting  points,  U0  =  5. 


Particle  Paths,  UQ  =  10,  VQ  =  -1 


Figure  9  -  Ion  trajectories  originating  from  four  different  starting  points,  U0  =  10. 
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Particle  Paths,  U0  =  20,  VQ  =  -1 


Figure  10  -  Ion  trajectories  originating  from  four  different  starting  points,  U0  =  20. 


Table  1  -  Initial  and  final  values  for  ion  trajectories  in  Figure  8  (U0  =5). 


# 

Co 

Po 

U0 

Vo 

00  O 

Up 

Vp 

0f(°) 

1 

-20.2 

3.28 

5 

-1 

-11.3 

5.10 

-2.94 

-30.0 

2 

-15.2 

3.66 

5 

-1 

-11.3 

5.13 

-2.79 

-28.5 

3 

-10.1 

4.25 

5 

-1 

-11.3 

5.12 

-2.07 

-22.0 

4 

-5.05 

5.26 

5 

-1 

-11.3 

5.08 

-1.20 

-13.3 

Table  2  - 

Initial  and  final  values  for  ion  trajectories  in 

Figure  9  ( U0  =10) 

# 

Co 

Po 

U0 

Vo 

0o  n 

Up 

Vp 

0f(°) 

1 

-20.2 

3.28 

10 

-1 

-5.71 

10.10 

-2.78 

-15.4 

2 

-15.2 

3.66 

10 

-1 

-5.71 

10.11 

-1.97 

-11.0 

3 

-10.1 

4.25 

10 

-1 

-5.71 

10.08 

-1.39 

-7.8 

4 

-5.05 

5.26 

10 

-1 

-5.71 

10.04 

-1.09 

-6.2 

Table  3  -  Initial  and  final  values  for  ion  trajectories  in  Figure  10  ( U0  =20). 


# 

Co 

Po 

U0 

V0 

00  O 

Up 

Vp 

0f(°) 

1 

-20.2 

3.28 

20 

-1 

-2.86 

20.07 

-1.69 

-4.82 

2 

-15.2 

3.66 

20 

-1 

-2.86 

20.06 

-1.37 

-3.90 

3 

-10.1 

4.25 

20 

-1 

-2.86 

20.03 

-1.17 

-3.35 

4 

-5.05 

5.26 

20 

-1 

-2.86 

20.01 

-1.04 

-2.98 

IV.  Conclusion 

By  applying  a  2D  sheath  solver  to  a  region  near  a  dielectric  corner,  an  estimate  of  sheath  thickness  and  potential 
could  be  obtained.  Using  the  calculated  potential,  one  could  predict  the  change  in  angle  of  an  ion  trajectory  due  to 
the  radial  acceleration  provided  by  the  sheath.  Depending  on  the  ion’s  initial  axial  velocity  and  position,  a  change  in 
angle  of  up  to  10°  could  be  seen,  which  is  not  insignificant.  Although  the  changes  found  by  the  model  were  not 
large  enough  to  explain  experimental  results  showing  populations  of  high  velocity,  high  angle  ions  in  Hall  thrusters, 
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the  data  suggest  that  the  sheath  could  still  have  a  noticeable  impact  on  the  way  ion  trajectories  evolve  in  the  thruster 
plume. 

Since  the  model  suggests  that  the  sheath  has  the  ability  to  influence  the  behavior  of  the  plume,  further  work  is 
warranted.  First,  as  discussed  in  Section  II,  the  2D  solver  should  be  developed  to  the  point  where  it  can  be  applied 
to  the  entire  region  near  the  dielectric  corner,  not  just  in  the  “quasi”  2D  fashion  described  in  this  paper.  This  would 
require  either  a  variable-width  solution  domain  or  a  complete  solution  which  represents  both  the  pre-sheath  and 
sheath,  but  would  lead  to  a  more  accurate  estimate  of  the  sheath  potential  profile.  Additionally,  a  wider  range  of  ion 
trajectories,  corresponding  to  an  array  of  thruster  operating  conditions,  should  be  studied.  It  would  be  useful,  for 
example,  to  see  what  effect  the  sheath  has  if  the  thruster  is  running  at  off-nominal  versus  nominal  conditions. 
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